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NOMENCLATURE 


Equation  ( 3 . 37a , b ) 


=(l/2Kot/RJ(t))-'/>,  constant 


Duany  variable  for  integration,  interface  spocd 
Temperature,  temperature  vector 


Velocity  of  the 
Thermal  diffusivity 


Source  heat  release  t 


j Source  point  location 


conditions  that  might  lead  to  an  exact  solution  of  the  Stefan 
interface  position  at  small  time  can  be  determined  graphically. 


position  of  time.  The  method  has  also  been  applied  for  testing 
element  method  has  been  used  to  verify  the  results  p 


distributions,  and  hysteresis  of  energy  storage  and  release, 
among  others,  which  have  not  been  previously  reported  in  the 


CHAPTER  I 
INTRODUCTION 


complexity.  Basically  a heat  and  masi 
that  is  prescribed  for  the  medium. 


e regions  where 


Id  tlie  solution  of  the  phase-change  problem, 

the  interface  (Stefan)  conditions  are  not.  Th« 
position  must  be  determined  by  using  these  nonlinear 

solutions.  Only  a feu  attempts  have  been  successful,  and 
problems  that  can  be  solved  exactly  have  been  limited  to 
that  are  adaptable  to  a similarity  transformation  (1-9]. 
these  problems,  the  interface  position  varies  as  a square  roi 


phase-change  literature  heavily  populated  with  approximate 
solutions.  In  these  efforts,  it  is  difficult  to  develop  an 
analytical  nethod  that  is  uniformly  accurate  over  a vide  range  of 
time;  see  Chapter  II,  Revicv  of  Literature,  for  a comparison  of 


resolving  the  singularity  at  the  phase-change 
solving  the  problem  in  a time  varying  domain, 
tend  to  incorporate  analytical  and  numerical 
a single  analysis  for  improving  the  accuracy  of 


y Lightfoot  [10]  f 


the  solution  of  the  phase-change  problems.  The  problems 

effects  being  negligible.  In  the  first 
i method  (SSH)  is  developed.  This  method 

or  multiple-phase  melting  and  solidification  problems  in  a 


In  the  second  method,  a boundary  element  method  (BEM)  is 
developed.  This  method  is  tested  only  insofar  as  it  is  related 


■ rarely  been  tested  in  the  literature  but  are 
•rtant  in  energy  storage  and  material  processing,  among 


provide  the  theoretical  basis  for 


the  physics  of  the  phase  change, 
gained  in  this  analysis  will  be 


s expected  that  the  insight 


exactly  remain  those  that  solved  by  Neumann  over  a century  of 


2.2  Literature  on  Exact  Solutions 
len  veil  established  today  that  for  a phase-change 

to  a similarity  transformation  [1-9].  This  limits  the  problem  to 

properties  and  imposed  with  constant  temperature  boundary 
conditions.  For  the  convenience  of  classification  of  problems 
reviewed  in  this  work , those  that  can  be  solved  exactly  will  be 

Stefan  problems  or  moving  boundary  problems,  following  the 


functions  [12-15]. 
as  a linear  relation  [12]  a 


[Id]- 


also  by  using 
9 modelled  as  a single 


Boley  [17]  introduced  another  unknown  for  the  heat  flux  at  the 
r problems  [18,19].  Boley’s 


problems  in  multiple  dimensions  [20]. 


indeed  the  case  with  the  Stefan  problems.  In  the  solution  of  t' 
Stefan  problems  by  the  asymptotic  expansions,  two  limiting  can 
are  derived;  one  dropping  the  transient  temperature  terms  in  t 


transient  interface-position  terms  in  the  governing  and  interface 

valid  at  small  time.  The  former  is  also  known  as  the  quasi- 
steady solution,  and  the  latter  as  the  quasi-stationary  solution. 


n problems  [21-28]. 


singularity  associated  with  phase  degeneration.  They  can  also  be 
extended  for  the  solution  of  Stefan  problems  where  convection  and 


method  of  weighted  residuals,  the  integral  acthod  follows  the 
method  of  von  Korean  and  Pohlhausen  that  was  originally  developed 
for  boundary  layer  analysis  in  fluid  mechanics.  Goodman  [3]  used 
the  integral  method  to  study  heat  transfer  that  includes  the 
solution  of  ablation  and  Stefan  problems.  Bell  [29]  refined 
Goodman’s  method  by  dividing  the  temperature  variable  into  equal 
intervals  and  derived  a system  of  first  order,  nonlinear, 
differential  equations  for  a set  of  penetration  variables  which 
defined  the  positions  of  the  isotherms  created  by  the  division. 
Charach  and  Zoglin  [27]  incorporated  the  integral  method  with 
perturbation  theory  to  study  the  solidification  of  a slab  that 
was  initially  superheated  but  had  a small  Stefan  number.  Vao  and 
Prusa  [9]  improved  Goodman's  results  by  incorporating  coordinate 
transformation  in  integral  analysis.  They  demonstrated  that,  for 
a material  with  a Stefan  number  less  than  unity,  the  integral 

than  the  method  without  transformation.  The  value  of  the 

known  as  classical  formulation.  In  the  classical  formulation  one 

consist  of  a distinct  interface  position  together  with  se 
functions  representing  temperatures  in  different  phases  a 
medium.  A formulation  in  marked  departure  f 


functions  are  used  for  the  enthalpies  in  different  regions  of  the 
■ediuo.  These  enthalpies  satisfy  an  integral  equation  derived  by 
using  the  original  beat  diffusion  equation  multiplied  by  a test 
function  and  integrated  throughout  the  time  and  space  domain  of 
analysis.  The  test  function  must  bo  selected  so  that  it  decays 
rapidly  as  the  space  and  time  approach  infinity — a condition 


The  principal  advantages  of  the  weak  formulation  are  the 
interface  position  vbich  can  be  determined  after  the  solution  is 
obtained.  In  addition,  in  the  weak  formulation  the  interface 
must  not  be  a sharp  boundary;  a mushy  extended  zone  may  also  be 
analyzed.  Most  of  all,  the  method  is  attractive  in  solving 


[30-32]. 


4 Approximate  Numerical  Solutions 


a stretching  and  translation, 
sing  an  analytical  expression 


translated.  Using  coordinate  transformation  leads 


the  further  nse  of  Fourier  modulus  for  the  normalization  of  time, 
the  original  Stefan  problem  can  be  transformed  to  a dimensionless 
form  that  can  be  solved  in  a fixed  domain.  Vievcd  from  this 
perspective,  the  method  is  akin  to  the  embedding  method  reviewed 
earlier.  Yet  normalizing  space  and  time  individually  permits  the 


Coordinate  transformation  itself  does  not  solve  a Stefan 
problem.  It  only  provides  the  convenience  of  working  in  a fixed 
domain.  Its  use  in  the  integral  method  to  improve  accuracy  was 
discussed  earlier,  and  numerous  studies  provide  examples  of  use 
in  numerical  work.  For  example,  Yang  [33]  used  it  in  the 

generalized  the  one-dimensional  form  of  the  transformation  to  the 
solution  of  a two-dimensional,  diffusion  controlled,  phase  change 
problem  in  a cylindrical  geometry.  Their  work  was  further 
extended  by  Saitoh  [35]  in  consideration  of  the  shape  of  the 
container.  Sparrow  et  al.  [36]  applied  the  transformation  to 
determine  the  effects  of  subcooling  on  a cylindrical  geometry  and 
verified  the  results  of  earlier  work  by  Tien  and  Churchill  [37]. 

transformation  in  a spherical  geometry  was  studied  by 
and  Kucera  [38].  Furzeland  [30]  mn< 


tuo  using  fixed  grids;  the  Utter  wore  found  to  produce  inferior 

that  vas  more  accurate  than  the  fixed  grid  aethods  vas  developed 
by  Murray  and  Landis  [39]. 


2.5  Recent  Solution  Methods 
even  today.  Some  recent  efforts 


The  Stefan  problems  have  recently  been  solved  as  a system  of 
region  [40]  and  finite  slab  [41]  subject  to  temperature 


if  order  approx i nation  o 


A conversion  of  the  Stefan  problem  to  an  infinite  number  of 
ordinary  differential  equations  in  tine  vas  also  attempted  [42]. 


phases  for  spatial 


and  Greif  [43].  In  this  approach, 
phases  vas  specified  as  a known  singl< 
analysis  vas  then  applied  f 


Frederick  and 


ngic-phase  solution;  an 

generation,  variable  properties,  a 
nonplanar  and  multidimensional  systems. 

Recent  analytical  studies  concentrate  on  the  exteus 
generalization  of  previous  efforts.  For  example, 

Greif 's  analysis  cannot  be  applied  t 

[25]  vas  presented  by  Gutman  [44].  Similarly,  • 

conductivity  and  specific  heat  that  vas  b 

Sunderland  [46].  An  analysis  using  s 
developed  by  Prudhomne  ct  al.  [47]  for  the  solution  of  problems 
in  plane,  cylindrical,  and  spherical  coordinates  with  the 
boundaries  imposed  vith  three  types  of  boundary  conditions.  The 


dependent  thermal 


Stefan  problems 


surface  flux  conditions  that  would  produce  a prescribed  interface 
motion.  in  this  effort,  the  problem  was  treated  as  an  inverse 
problem  [48]. 


efforts  arc  focused  on  the  use  of  boundary 
Primarily  an  integral  equation  set hod,  the 


[49]  and  the  temperature  field  associated  with  freezing  around  a 
buried  pipe  [50].  A two-diaensional  solidification  problen  with 
an  iterative  implicit  algorithm  to  determine  the  interface  motion 
was  developed  by  Zabaras  and  Uukherjee  [51].  An  elementary, 
noniterative  system  featuring  linear  interpolation  over  elements 
on  a polygonal  boundary  was  developed  for  the  solution  of  two- 
dimensional  problems  by  O'Neill  [52].  The  boundary  element 


Stefan  problems  in  this  dissertation.  The  present  effort  will  be 
focused  on  the  evaluation  of  the  accuracy  in  these  methods. 


monographs  and  texts,  which  will  be  mentioned  here  1 

Yeo  and  Prusa  [9]  reviewed  tbe  solution  of 
freezing  problem,  in  a recent  issue  of  tbe  “Advance,  in  Heat 
Transfer."  In  this  monograph,  literature  on  Stefan  problems 
dominated  by  convection  and  radiation  effects  are  reviewed 


mograph.  Literature  c 


with  tbe  classical  single-phase  problems 
mathematical  difficulties  inherent  in  I 


entitled  “One-dimensional 

Crank  [7]  published  a book  entitled  “Free  and  Moving 

the  solution  of  both  Stefan  problems  and  the  allied  free  boundary 
problems.  The  book  provides  a broad  perspective  into  the 
solution  of  problems  in  a time-variant  domain.  While  the  Hill 
book  only  covers  the  analytical  solutions  to  Stefan  problems,  the 
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phase,  imposed  with  cyclic  conditions,  the  solution  is  dividt 


W). 


To(*,0)=0  (3.2) 

To(0,t)=F(t)  or  -h'>r''^,t)  =G(t)  (3.3a,b) 

T„(oo,t)=0  (3.4) 


uatioo  (3.2)),  t 
tine  aero,  0<F(0)<T, 


nelting  and  0>F(0)>T„  for 


Second  Stage 


Tt(0,t)  = F(t)  or 


a(*!*o)  — Fo(a,ta) 


(3.8) 


T.(R,(t),t)=Tm  = VR1(t),t) 


W.(R^t).t)  OTtdMt),t)_^L  dR,W 

R,(to)=0 


^=i  ^ . Tc(x,t), 


(3.13) 


3!T„  , 8Ta 

*?"~s  "Jr  • 


T.(x.t), 


T.(x.t)  =T„(x,t,)  for 
Boundary  condition: 


*,(t)<*<00 


(3.18) 


=T„=Te(R,(t),t) 

”,lW,t)  art(R,(t),t)  Pt  diyt) 


(3.20) 

(3.21) 


T„(R,(t)  ,t)  =T„  = VR,(t),t) 

*T,(M t) .t)  OTt(R^(t),t)  pL  dR,(t) 


changed  if  HOP  precedes  HIP. 


sides  of  (3.21) 


[10].  For  tho  | 


*l(to) : 


R,(t.)  = i 


(3.31) 


Initial  conditions: 


T(*,t)=T0(x,t„)  for  t = tg 

T(x,t)  =Tj(x,t|)  for  t = t, 

T(x,t)=T,(x,te)  for  t = te 


(3.32) 

(3.33) 


(3.36) 


deleted. 


n solving  the  second  s 


in  (3.25)  a, 


equations  (3.29),  (3.31),  (3.32),  (3.3d)  and  (3.35)  must  be 
deleted.  Finally,  when  solving  the  third  stage  problem, 
equations  (3.30),  (3.32),  and  (3.33)  most  be  deleted.  Again,  the 
sign,  preceding  «,  and  must  be  interchanged  if  HOP  precede.  HIP. 
For  the  temperature  cycle  imposed  on  a subcoolcd  solid  as 

will  eventually  merge  to  the  effect  that  the  entire  medium 


roiers  to  the  time  of  merging  two  fronts.  After  te,  A,  and 
term,  in  (3.26)  disappear,  and  cooling  in  the  solid  follow,  t 
conditions  given  by  (3.26),  (3.27)  and  (3.36).  It  is  noted  tha 
for  a solid  that  is  not  subcooled  initially,  this  fourth  staj 
will  not  occur  as  will  be  shown  later  with  examples  for  bo: 
cases  in  this  dissertation. 


2 Solution  of  First-Stage  Temperature 


’ > »=°  (t-s)1' 

fJL  Iif2 
( ® *-  1 
l £ C(s) 


solving  the  following  two  equations: 

F(»o)=T„ 

ic  temperature  is  imposed  at  the  surface,  and 


when  the  heat  flux  is  imposed  at  the  surface.  Clearly,  this 
equation  (3.36). 

Two  special  cases  are  worthy  of  note.  Vhen  F(0)>T„jS0  for 
a melting  problem  or  F(0)  < T,n  ^ 0 for  a freezing  problem,  phase 
changes  immediately  upon  imposition  of  the  boundary  conditions. 


still  two  phases  in  the  medium,  and  the  problem  is  known  as  a 
two-phase  problem.  Another  case  arises  when  T„(x,0)  = T„  = 0, 
corresponding  to  a medium  without  superheating  or  sobcooling. 


i:.. 


to)  k J^=to  C(x,t|K,(r).r)  d: 

— *«»-{;  - s* 


x = B,(t)  and  T,(R,(t),t)=T„ 


(3.45) 


T„  = T0(R,(t),t) 

*»(*-«.)  if  dr 


T,(x,t,)=T0(x,t,) 

* fe  ^ C<*.‘.l«,M.r)  dr 

(3.46) 


T,(0.t,)  = F(t1)=T„  . 


T„  = T„(0,t,)d:  \l  0(0,t,|»1(r),r)  dr  , t,  > t„ 


I~=o  T,(x,,t,)  G(’‘-t‘t'|x'’0) 


+ k J~‘  [-dR|J?t,)  ®(x»t*t||RI(r+t(),r) 


~ J~=oJ«=«o  m 

= - J^=to  C(K,t|R1(^),^>dr 


Ta(X,t)=T0(x,t)+fe  (-^)  C(x,t|R|(r),T)  dr 


G(x,t|R,(r),r)  dr 

(3.82) 


,=R,(t)  and  T, in  (3.52)  a. 

T.  = T„(R1(t),t)+i  P=tjj 

+jf  J^=t,  ediW.tlM*-).’-)  * 


Tm  = T0(Ra(t>,t)+fe  (T^^)  C(R,(t),t|R1(r),r)  dr 

«(R,(t),t|RJ(r),r)  dr 
(3.54) 


T,(x,t)=T0(*,t)+fe  ( 

*s  j:u, 


integrals 


effective.  Thun  for  a cyclic  condition  imposed  on  a medium 
without  initial  subcooling  or  superheating,  the  interface  R, 

duration  of  the  cycle.  As  a result,  the  relevant  integral  for  R, 
(first  integral  in  (3.52))  vould  run  from  r=0  to  rat,  whereas  the 
relevant  integral  for  R,  (second  integral  in  (3.52))  would  run 
from  r=t,  to  r=t.  On  the  other  hand  for  a monotonic  condition 


t and  (3.52)  r 


Nothing  compromised  in  the  derivation  above,  the 

expressed  in  an  integrodifferential  form.  A numerical  method 

time  being,  equation  (3.43)  will  be  tested  for  some  boundary 
conditions  that  might  lead  to  an  exact  solution  of  the  problem. 


Use  is  eade  of  the  following  for  the  boundary  temperature 


F(t)=(C,+C,t)>T„,  R,(t)=C3t,'!  (3.56n,b) 

Substituting  then  into  (3.45)  lends  to  three  relations  as 

*o  = ® (3‘87) 

erfc^sgl;  p exp(Pa)erf  (P)erfc(P)] 

" (3.58) 

2 ^ [(J+l»)erfe(P)-£  exp(-PJ)]=0  (3.59) 

where  P represents  one  half  of  the  reciprocal  of  the  square  root 
of  the  Fourier  modulus  defined  ns  ot/H,,(t),  and  Ste  represents 
the  Stefan  number  which  is  defined  as  eT^/L. 

To  satisfy  (3.59)  at  all  time,  either  C,  or  the  bracketed 
expression  must  be  zero.  For  that  bracket  to  vanish. 


exp(PJ)  J“exp(-y*)  dy  = 


(3.60) 


(P+^Ps+2)  - 1 < exp(P!)  f”exp(-yJ)  dy  < (P  + 'JP’tj)  1 

(3.61) 

for  P>0,  according  to  Aramowitz  and  Stcgun  [83).  Use  of  (3.60) 


P’  + P^4<2P’  + .<I*  + ^ 


G(t)=C4+Cst,  R,(t)  = C6t,  T„  = 0 


C^+Ar  C,)t+(^+^,c,s)t,+  - • ■ - =o 
<£  + & C«)t,/,  + (w  + ^ + T2S5  C*’)t3/,+ 


c,=  -!£c„,  CS=-J^C„3 


It  follows  that  the  coefficients  C , and 
of  C„  and  C,s,  respectively.  Now  since  I 
advancing  melting  front,  C,  and  Cs  mus 
(3.67),  a result  in  violation  of  ’ 
that  is  imposed  at  the  boundary. 


ust  have  oppsite  signs 
negative  according  to 
rulos  out  the  possibility 


interface  position  i 


and  interface  position  combinations.  Since 

of  the  right  hand  side  of  equtaion  (3.45)  while  the  interface- 
position  effect  has  been  confined  to  the  second  term  of  the  right 
hand  side  of  this  same  equation,  the  analysis  above  can  be 
extended  for  testing  new  combinations,  such  as  (3.56a)  and 
(3.64b),  and  (3.64a)  and  (3.56b).  It  can  bo  shown  that  the  exact 
solutions  are  impossible  for  these  combinations.  The  tests  given 
above  show  clearly  that,  other  than  a constant  temperature 
imposed  on  a semi-infinite  medium  with  uniform  initial 
ot  and  linear  interface  relations  arc 


8 Numerical  Solution  of  Interface  Positions 


be  solved  numerically.  The  numerical  solution  is  effected  b; 
moving  the  dR,(r)/dr  in  (3.45)  and  the  dH,(r)/dr  and  dRj(r)/dr  ii 
(3.53)  and  (3.54)  out  of  their  integrals.  The  interface  positio: 


±S(tNl-t,)  eCH,(tN|),tNll»,(r),r)dr) 

= Jj  [T„-T0(R,(tKl),tNi)]  (3-®S) 

R,(t»)=0;  ",=1,2 N,i  tWj<V 


= £ [T„-T„(R,(tNj),tKj)] 


In  a similar  fashion,  (3.54)  can  be  written  as 


dr) 


.^(^(t^J.t^l^CrJ.Odr) 


=C  tT„-T0(R,(tWj),tWi)] 


(3.70) 


where  R}(tu  )n0;  N,=N,+1,  N,+2,...,Sjj  t,  < tN  <t,,  or  t,  < ts  <tc. 
Equations  (3.69)  and  (3.70)  will  be  solved  simultaneously  for 

gradient  dR/dt  arc  related  as  shown  in  Figure  3.2.  The  position 

(3.43),  (3.52),  and  (3.65)  to  find  T(x,t).  Finally,  the 

the  following  four  equations: 

Jo"  q(t)dt  = pj“  cT(x,t„)  dx  (3.71) 


q(t)dt  = , {|*l(t")[l-  + eT(x,t„)]  dx 


tor  toSt.st,, 


Jo"  9(t)dt  = p (JJl(t")[L  + cT(x,t„))dx+  J“(tn) 
-J^Ldx) 


q(t)dt  s , 


cT(«.t„) 


(3.74) 


for  te<t„<tp.  They  are  derived  strictly 

heat  from  the  sensible  heat,  an  important  par' 
Stefan  problem.  The  solution  is  now  complete. 


An  SS*  computer  program  useful  to  solve  a Stefan  problem  i 
a subcooled  medium  imposed  with  a cyclic  temperature  condition  a 
that  shown  in  Figure  3.1  is  provided  in  Appendix  A. 


consideration  is  given  to  the  fact 
carried  out  in  regions  whose  doi 
time  (Figure  4.1),  and  (ii)  the  c 


boundary  element  methods, 


of  these  domains  are  unknown  a priori  and  mus 
part  of  the  solution.  These  characteristics  distinguish  the 
formulation  developed  in  this  chapter  from  those  of  the  solution 
of  conventional  heat  diffusion  problems  in  which  the  domains  and 

4.1  Derivation  of  Boundary  Integral  Equation 
Boundary  integral  equation  for  the  solution  of  a moving 
boundary  problem  will  be  derived  on  the  basis  of  solving  a two 
dimensional  beat  conduction  problem.  There  is  no  heat  source  in 
the  medium,  the  material  properties  being  constant.  As  shown  in 

change  interface.  For  each  phase  indicated  in  the  figure,  the 
s denoted  by  0;  its  boundary  is  T,  both  being  functions 
The  governing  equation  in  each  phase  is  given  by  the 


Figure  4.1  Phase  regions  used  in  the  derivation  of  the 
integral  equation  in  BEM  analysis 


HU  *+<  IIW>  ’■*  S “•  «■*> 

Then  by  using  the  Green’s  second  identity,  the  left  hand  side  of 
equation  (4.2)  can  be  written  as 

Ilia.) « “-IIIto  I'S-’ffi  "* 

+HJa.)TV’r  dn,lc  <*•» 

shown  that  (52) 

JiKt)  &(T‘T)  dn  = &In(t)  T'T  dn- Jnt)  T*T(S'7)  dr  (4  4) 

where  1 denotes  the  outward  drawn  normal  to  the  boundary,  and  v 
With  the  use  of  (4.4),  one  derives  an  alternative  expression 

Ilia.)  rl  U JUa.)  1 “» 

= " lau)  1 rT  Ilia.)  * ”))  d!wc 

-JllmiS’-'i1'5)'™  «■•) 


it  vanishes  at  time  t.  Using  equations  (4.3)  and  (4.5)  in  (4.2) 


(4.6) 


R'ln.)  drdl+"  law n 

t can  be  shown  that  T*  is  the  solution  of  the  following 


oV’T*((,tF!a,t)+^(£.tF!X,t)=  -£(£-x)<(tp-t)  (4.7) 


which  represents  the 


functions.  Using  (4.7)  in  (4.6)  gives 


’«.v)-J^(rW  [rg-'C*lrr(r.vi] 


(4.8) 


ad.  the  point 


«!»«.«— (JJn.)  drd‘ 

*Jwr,do  (...) 

where  C(£)  depends  on  the  location  of  the  point  £ as  well  as  the 


boundary  [54-66]; 


S-»S- 


T„(x.t0)=T0(*,t0) 


(4.15) 


to  reach  a phase-change  temp 
discussed  in  Chapter  III. 

Application  of  (4.9)  to 


0«)T(«.tp)=“|t^  [T,((.tFin,.t)OT(^’t^ 

-T(H,,t)”t«,yl>»-t>^  T(R,,t)T*({,tF;R,,t)^ 

_T.(£itp!0,t)2I^+T(0,t)«C«^i^]  dt 

(4.16) 

the  boundary  0 is  stationary  and  %(tg)  =0.  Similarly, 


C(£)T«,tF: 


+ J0  T»«,»F»x,1*)T,(*,t0)  dx 

ements  over  vhich  the  temperature,  flux  and  interface 


*s  ft, 

<*■»> 


- jjt, 


in  the  integrands.  Using  (4.20)  in  (4.18)  and  (4.19)  leads  to 


Cjf  = 4 jztff  i '<>’  T»+i  <^)/  T- 

+i  fett  4*fgu  n • ‘“O-1  (4 


cflT' =4 <4-4^  T--i  1,“"' (^5'  T" 

♦IT* 


H&'  = H,y  + PA/,F  • 

n> 


i={  ®-*  i.  Eg;  £ «•»> 

tL  «•*> 


= -C0,  P,=C,  for  region  b,  nod  P,  = 


(4.27) 


Region  b~ 

Tf  = -£  J;Cg'  Mfi-fanr  T.+i  <^)'  T„ 


Region  •- 

xf  = J q{..+J:fiS'  T„-i  C^>' 


Expressions  for  T* , ffT*/dx,  and  C,*  can  be  found  in  texts 
f heat  diffusion  problens  [54-56).  The  Green’s  function  T*  in  a 


^Ftj] 


(4.31) 


that  arc  positioned  at  cither  side  of  the  point.  Referring  to 
the  two  dimensional  domain  shown  in  Figure  4.3,  it  can  be  derived 
a.  [55] 

Cf  = £ + = <«-32) 


to  be  taken  as  unity  for  the  sense  point  located  interior  to  a 


reaches  the  phase-change  temperature.  This  temperature  has  been 


1-  i^k  C.  ■—[-Kivs 


I"o  * '■ 

■•^-EifeT,-tl-*'lJWvS^i>I  ‘*'W 


e similarity  of  one 


term  in  equations  (4.36)  and  (4.37), 
facilitating  the  development  of 


integration  oust  still  be 


hen  for  an  arbitrary  t. 
. be  R.Mmpt  + q,  wl 


”*(°ft,0’t)  dt  = 0 (4.38.) 

OT,(°^i0,t)  dt  = ° t4'38") 

= , W‘(0'^'il>l’t)  dt=  cnp[p(pyi>] 

^__JjP(tr-t/)*(ptr+*i) j _ ■■■.tP(tF-t/-i)*(Ptf-t|i))  J 


(4.39a) 


, «T-(0.^,R,Ja  „„  .,,[^31, 

[i  -ft,(t-t-i);(ptr^]  (4'39b) 


-i-i— "<j,.K-gb1]  <*■“> 

•si-w,  - 

= exp[_P(Mi£2^£l2iI].|“j_i  £ „„(-«»)  du 

(4.41a) 


It  = ^-crf[— 


«^=0f^_,  T*(o,tF!o,t)  dt=jf 

(4.4! 

C«F  = “Jt£_,  T*(0,tF!0,t)  dt  = Jf  gtp-tF., 

(4.41 

Gm,  = 0{t{_,  T*(0,tf.iH,,t)  dt 


(4.43a) 


= 5l<j£>  [e*pC-«J., -«})/“/ 

- 'ft*(erf(uy)  -erf  (uy_j)))  (*.«*, 


= “JtF-I  T'(R"tfi°’t)  dt 

[«xP<-4_1)/ur.1- 


55  L 1 J4o(tf-ty)  Jte(VV-i) 


(4.45b) 


interface  position 


stationary, 


Equations  (4.1 

«£'=» 

N (4.46a) 

8£r=‘ 

fi £'=, 

Ri(tp)-Ri(t|) 

" 1 ' [Crf(  j4o(tF-~>  ",<j4o(tF-t— ' J 

(4.47a) 

Bf.F  = 

|tp  «T"(R,,tFiR,,t)  j. f4.47bi 

G„V  = 

[exp(-v*_,)/v^_I - esp(-vj)/vy 
-fT- (erf(vy)— erf(vy-i))] 

R,(t,)  j Ei(t/> 

= Rl(V)^,Rl—  [«XP(-C 
-exp  (-(»/-v,)J)/(«/-v/) 

-q»(erf(»,-v, )-«*(:!,_, -V/_i)>] 


"/-'■jtoCtf-t,.,) 


(4.49a) 


, the  following  equations  are  derived: 


s fig/  by  changing  R|(tj?)  to  (( ty)  (4.60a) 


as  «"■  by  changing  R,(tp)  to  f(tF)  (4.50b) 


r[  P«(tr)-jPtr>l))]  J°/_i  £ exp(  - u3)  du 

(4.51a) 


■■is., 1 


.[-I— 


of/  by  changing  R,(tF)  to  «tF)  (4-52«) 


f by  changing  Ri(tF)  to  f(tF) 


,g( tF)-Ri(t/-i).  _2gf'/1 
j4o(tF-t,_,)  11  ^ 


Again  for  a stationary  interface. 


Og/  = ojij _t  T*({,tF!R,,t)  dt 

= {(V)^Bl(^)  [cxp(.(y/.V/)t)/(y/.1-v/.1) 
-expf-Cyy-vyj’J/fyy-v,) 

- ■ (erf  (y,  - v,)  ) - erf  (y,  . , - »,  -,))] 


(4.84a) 


=«V)-R.I^  [exp(.(yF., 


-irt/typ-t-XF-i) 


4.7  Numerical  Solution 

A hybrid  method  L used  to  .olve  the  Stefan  problem  by  the 
boundary  element  method..  For  a nubcooled  or  nuperheated  medium, 
the  temperature,  in  the  first  .tag.  are  .olved  analytically  a.  in 
Chapter  HI.  These  analytical  solutions  are  then  used  as  input 
to  evaluate  the  T&„  as  shown  in  Section  4.4. 

Next  equations  (4.28)  and  (4.26)  are  structured  in  a matrix 


C.Q-E-T=0 

where  6 and  H represent  coefficient 
elements  G and  H whose  expressions  have 


n derived  previously. 


The  q and  T denote  vectors  of  heat  flux  and  temperature, 
respectively.  For  example,  Q contains  q{,  qf,t,  and  q{,„,  vhilst  T 
contains  T T{>  = T„,  and  Tf0  = Tm,  where  f=l,2,  --,F.  On  the 


right  hand  side,  £ denotes 


Equations  (4.25)  and  (4.26),  as  derived,  signify  t 
e integration  over  time  starts  from  tq  up  to  tp,  the  sent 
ius  the  numerical  solution  involved  is  still  a time  l 

.quires  the  us.  of  those  at  time  F-l  and  earlier  for  inp 


calculated  for  each  sense-* 


II  and  only  one  nev  G must  be 
oint  pair,  a property  greatly 


n numerical  computation. 


temperatures  at  the  interior  points.  In  addition,  calculation  of 


method  particularly 


advantages  collectively 


4.8  Computer  Program 

A BEM  computer  program  useful  to  solve  a Stefan  problem  in  a 
semi-infinite  medium  with  or  without  subcooling  or  superheating 
and  imposed  with  constant  or  monotonic  temperature  or  flux 
conditions  is  provided  in  Appendix  B. 


DEVELPOMENT  OF  CHARTS 
CONVERGENCE 


CHAPTER  V 

ND  TESTS  FOR  UNIQUENESS,  ACCURACY, 
iD  STABILITY  OF  SOLUTIONS 


The  structure  of  equation  (3.68)  is  worthy  of  note  it  is 
well  suited  for  the  development  of  computer  algorithms.  Clearly, 
the  right  hand  side  of  this  equation  (second  term)  hinges  on  the 


right  hand  side  of  this 


intcgrablc  in  closed  form  [53,57].  These  characteristics  enable 
the  algorithms  to  be  developed  for  an  accurate  solution  of  the 
Stefan  problem. 

To  folly  utilize  the  analysis  developed  in  this  work,  some 
charts  are  developed  for  the  solution  of  the  interface  positions 


Chart  for  Evaluation  of  Interface  Position  for  Boundary 

For  a semi-infinite  medium  of  constant  properties  imposed 
with  a constant  temperature  condition,  the  interface  position 


Stefan-Neumann  problem,  for  which  the  interface  motion  can  be 
obtained  by  using  (3.58)  which  is  plotted  in  Figure  5.1. 

temperature  imposed  on  the  surface  and  the  Stefan  number  is 


contain  the  surface  temperature,  which  has  been  grouped  together 
with  the  phase-change  temperature  Tm  to  form  a temperature  number 


for  either  melting  or  freezing  according  to  the  coordinates  shown 
in  Figure  3.1.  For  a one  phase  problem,  Tm  would  be  zero.  Then 


position  must  generally  be  solved  numerically  as  discussed  in 

The  interface  position  at  small  tine  is  estimated  by  setting 
N[  to  unity  in  equation  (3.68)  as 


St.-cT„/L 


-gfegg  )=±[1-M*^i*i]  (5.1) 


interval,  Fol=e(At)1/R13(tI),  and  SSF  and  SSG  represent  the  braced 


on  the  boundary. 


see  SSF  and  SSG 


he  plotted 


P(t)=I  [®xp(S^)-l] 


Properties  for  Materials  Tested 


interface 


This  problem  ia  often  considered  as  physically  insignificant 


eqivalent  to  imposing  a 


criterion  is  taken  to  be  3xl0_l#,  and  the  results  (not  shoun) 


exponential  conditions  given  in  (5.2)  and  (5.3),  the  interface 
curvature.  As  shown  in  Figures  5.3  and  5.4,  the  interface 


Accuracy  test  of  the  SSM  in  the  solution  of  the 
interface  positions  of  tvo-phase  melting  of  aluminum 


close  to  0.001%  for  two-phase  melting  (Figure  5-4).  More 

The  same  tests  are  repeated  for  two-phase  melting  of 

R,(t)  can  be  calculated  and  plotted  in  Figure  5.5.  This  time  a 
d stability  of  the  SSM. 


example  (Stofan-Neumann  problem) 


suppressed 


one-  and  two-phase  melting  of  aluminum.  The  accuracy, 
convergence,  and  stability  of  the  BEM  are  thus  verified  as  well. 


vo-phase  Bolting  of  paraffin 


ERROR  OF  R,(t)  (X) 


of^two-phase  melting  of  aluainua 
t temperature  condition 


c Condition  of  Tine 


Trends  of  Interface  Position  Curves 

conditions;  see  Table  6.1.  The  results  are  plotted  in  Figure 
fl.l,  where  the  length  I,  used  in  the  dimensionless  groups  is  taken 


boundary  conditions 


e of  the  surface  temperature  d 


for  both  one-  and  two-phase  melting.  It  is  also  expected  that  a 
sensible  heating  prior  to  melting  delays  the  phase  change;  see 
curve.  7 through  9.  However,  the  trend,  of  these  curve.  re«,in 
e important,  the  interface  positions  for  curves  4 
y almost  linearly  with  time,  which  suggests 
s square  root  relationship  for 


that  an  indiscriminate  u 
the  interface  position  in  scaling 
coordinate  transformation. 

Discussion  of  the  trends  due  to  temperature 


■ subject  t 


5 now  directed  t' 


results  are  plotted  in  Figure  6.2.  Trends  of  these  curves 
summarized  as  follows: 

1.  By  examining  curves  1 through  3,  it  may  be  said  that 
flux  conditions  imposed  at  the  boundary  for  time  greater  t 
zero  have  an  effect  of  changing  the  slope  of  the  position  cur 


counterparts  of  the  curves  1 through  3 in  Figure  6.1.  However, 
large  step  input  of  heat  at  the  boundary  at  time  equal  to  zc 


!)• 

2.  * sensible  heating  prior  to  melting  delays  the  onset  of 
phase  change;  see  curves  5 through  7.  Again  the  slopes  of  these 
curves  are  consistent  with  the  rate  of  heat  input  at  the  surface. 
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Figure  6.2  Trends  of  interface  position  of  Stefan  problems 

imposed  with  constant  or  monotonic  flux  condition  of 


subcooling  h 


r over  a wide  range  of  time.  This  suggests  that  the 
as  an  effect  on  the  interface  position.  In  general 

c transformation  is  thus  justified, 
position  curves,  namely,  curves  6 and  8 in  Figure  6.1  and  curves 

the  literature  [9,11-14]. 

by  Goodman  (integral  method  (IK))  [3]  and  Ozisik  (Coupled 
integral  equations  (CIE))  [40].  It  is  found  that  using  a cubic 


= different 


CIE  resoles  are  good.  For  tl 
interface  position  is  found 

Figure  6.1.  Since  the  erro 


E method;  again  tl 


e methods  are  small,  they 


are  not  re-plotted  here. 

Vben  the  surface  temperature  changes  rapidly  w 

rapidly.  As  shown  in  Figure  6.3,  there  is  good  agreement  between 
the  present  method  and  the  CIE  at  small  time.  However  at  large 

as  much  as  20%,  and  such  error  is  expected  to  grow  further  with 
time.  This  completes  the  tests  of  the  temperature  conditions; 
attention  is  now  directed  toward  the  heat  ; 


functions  (CEF))  [12-14]  are  accu. 
slightly  larger  time  of  2 sc 


[11]  and  Tao  (complimentary  error 
ccurate  over  a small  tine  period  of 


m figure).  Nevertheless,  Goodman’s  IM  [3] 


.r  flux  condition  imposed  a 
2 and  Figure  6.5.  Again  E 


82 


0(0-3*10®  (W/m2) 


Figure  6.4  Comparison  between  the  SSM  with  the  acthodn  developed 


method  [11]  yields  results 


solution  method  at  small  time  (see  top  figure).  It  ag 
at  large  time  (sec  bottom  figure).  However,  Goodman's 

varies  rapidly  with  time.  Using  c 


I [3]  i 


be  found.  Since  the  strength 


Two  examples  are  used  for  testing  cyclic  conditions  imposed 

P(t)  =932  + 200  sin(^jt)  (K)  (6.1) 

G(t)  = 3«10»  sin(fljt)  (V/m>)  (6.2) 

where  the  period  is  20  sec ; see  Figure  6.7.  Only  one  cycle  is 
tested  for  both  examples.  As  before  (Section  5.3),  a bisection 


HUE  . I (•) 


in  the  literature  for  the  solution  of  Stefan  problems 


F(t)-932-200Bln(it) 

T(*,0)-932-^ 


1 1 the  overall  heat  input  over 


front  lags  slightly  behind  the  melt 
be  explained  by  referring  to  the 


Figures  6.9  and  6.10.  Figure  6.9  covers  times  from  1 through  5 
sec  and  11  through  15  sec.  These  times  cover  the  rapid 
temperature  rise  and  drop  parts  of  the  HIP  and  HOP;  see  Figure 
6.7.  Putting  them  in  one  figure  permits  comparison  of  the 


Interface  poaitiun  curves  for  a Stefan  problen 


Figure  6.10  Temperature  profiles  at  times  6 through. 10  sec  and 


parts  of  the  cycle.  This  leaves  the  rest  of  the  temperature 
curves  to  be  placed  in  Figure  6.10i  here  again  they  cover  the 

by  locating  the  point  of  intersection  of  the  curves  with  the  x 
axis  at  zero  temperature.  Then  according  to  Figure  6.9,  the 
freeze-front  positions  always  lie  slightly  to  the  left  of  the 
melt-front  positions,  indicating  that  the  freeze  front  lags 


front  at  10  sec  (Figure 


plays  in  delaying  the  f 


One  application  o 


e phase-change  material 
is  stored  in  different  phases  of  the  medium.  Use 

(3.72)  and  (3.73) i whereas  Figure 
the  difference  of  the  Q’s  at  t„ 

(*„-*•»-!)!  thc  r'sult  iB  then  “kc 


Referring  t 
ly  the  shape  of  the  la 


figur. 


is  plotted  by  computing 


he  advance  of  time,  resulting  in  i 
8 instantaneously  stored  (Figure  I 


se  of  the  sensible  heat  stored  i 


release  curve  follows  that  heat  store  curve  in  Figure  6.11  but  in 


As  mentioned  earlier,  this  gap  will  eventually  be  closed  moving 
Thc  sensible  heat  stored  lags  behind  the  surface  temperature 


(t)  (KJ/m2  ) 


Figure  6.  IX  Cumulative  heat  stored  or  released  for  a Stefan 

problem  imposed  with  cyclic  temperature  condition 


, (t)  (KW/m1) 


Instantaneous  heat  stored  or  released  for  a Stefan 
problem  imposed  with  cyclic  temperature  condition 


2 

t 


TIME  . t (sec) 


majority 


DISTANCE  FROM  SURFACE  . « (m) 


Figure  6.15 


Surface  temperature 


with  cyclic  flux  condition  without 


Instantaneous  heat  stored  or  released  for  a Stefan 
problem  imposed  with  cyclic  flux  condition  without 


only  one  cycle  is  tested  in  this  dissertation. 


equation  (3.24),  the  aaxi 


This  retracted  aclt 


ng  initially 


DISTANCE  FROM  SURFACE  . x (m) 


Temperature  profiles  at  6.7 


quantitatively  in  Figure  6.22  vhere  the  cu.ul.tive  heat  atorage 
and  release  are  plotted  for  the  phase-change  wax. 


are  plotted  with  curves  A,  B,  and  C,  respei 
Clearly,  curve  A is  generated  by  adding 
Referring  this  figure  vith  the  curves  for 
plotted  in  Figure  6.19  provides  insight 


stored  only  covers  the 


storage  is  quite  appreciable,  extending  all  the  way  from  time 
aero  to  t . There  is  a considerable  amount  of  heat  penetrating 


designers  working  in  the  field. 


For  the  temperature  cycle  imposed  at  the  surface,  the 


cumulative 


Q„(t)  (KJ/ml) 


Figure  6.22  Cumulative  heat  stored  or  released  for  a multi-phase 


the  sensible  heat  storage  continues  to  grow  after  this  point  in 


shown  in  Figure  6.23.  In  this  figure,  the  cumulative  heat  stored 

K (initial  temperature)  and  normalized  by  the  peak  temperature 

Figure  3.1).  This  plot  is  thus  analogous,  to  a certain  extent, 
to  the  hysteresis  curve  of  flux-density  versus  intensity  plot  in 

mainly  a manifestation  of  the  permeability  of  tho  magnetic 
heat  storage  and  release.  Clearly,  because  of  the  phase  change, 


scale).  But 


subcooled  medium' imposed  with 


energy  storage  and  release 


the  line  segment  OR  along  the 


CHAPTER  VII 

DISCUSSION  OF  RESULTS  BV  BEH 

Only  constant  or  nonotonically  varying  temperature  and  flux 
conditions  given  in  Tables  6.1  and  6.2  are  tested  by  the  boundary 
clcecnt  methods.  Again  aluminum  is  used  for  the  phase-change 
material  uhose  properties  arc  summarised  in  Table  5.1.  For  the 


expected  as  one  compares  Figures 


e digits  at  large  t 


o tvo  digits  at  small  tins 

Remaining  tests  vith  monotonic  temperature  and  flu 
conditions  yield  similar  trends  as  presented  in  Tables  7. 
through  7.6.  These  good  results  can  be  ascribed  to  using 


1 i terature 


[51].  Generally  in  the  BEM,  a starting  solution  is  necessary  for 

due  to  the  singularity  associated  with  the  condition  at  initial 
tine  [51].  Without  the  starting  solution,  the  BEM  nay  suffer 

variable  material  properties  can  be  easily  accounted  for  in  the 
output  (see  Section  4.7  in  Chapter  IV).  Whereas  in  the  SSM,  they 


uialyzed  in  Chapter  III.  For  such  problems, 
R3(t)>  and  Vet,  (4.25)  can  be  used  for 


method,  material  properties 


n the  solid  and  liquid  phases  have 


s boundary  element 


disappearance  of  the  interfaces  and  the  return  of  the  medium  to 

In  the  SSU,  a single  equivalent  problem  ((3.25)  through 
(3.35))  is  derived  in  which  a heat  sink  is  associated  with  a 
melting  front  while  a heat  source  is  associated  with  a freezing 


integrodlfferential  equations  is  derived  for  the  solution  of 
temperatures,  given  by  (3.36),  (3.43),  (3.S2),  and  (3.55) i and 


i an  accurate  determination  of  the  positions, 

condition  of  tine,  the  third-  and  fourth-stage 
solutions  do  not  exist.  On  the  other  hand,  for  a medium  without 


This  consists  of  weighting 
boundary;  see  (4.9).  Then 


moving  the  sense  point  to  the 


two  phases  in  the  medium;  see  (4.16)  and  (4.17).  In  the 

numerical  solution  by  the  BEM,  the  time  interval  is  discretized 
into  constant  elements  over  which  the  temperature,  flux  and 

derived  as  (4.27)  and  (4.28),  while  those  for  t£,,  (T,  and  G, 
(4.36)  through  (4.54b). 

imposed  on  a semi-infinte  medium  with  uniform  initial 
To  fully  utilize  the  SSM  solution  derived  in  this  work,  a 
imposed  with  a constant  temperature  condition.  Charts  for  the 


5.  Using  tuo  examples  that  include  a one-phase  melting  of 
tvo-phasc  melting  of  aluminum  and  paraffin  wax  imposed  with  a 
for  accuracy,  convergence  and  stability  of  solution.  It  is  shown 
positions  con  be  determined  with  good  accuracy.  For  example. 


The  SSM  are  thus  found  to  be  highly  accurate,  while  the  accuracy 


coordinate  transformatioi 


melting  front  continues  to  advance  e 


though  the  surface  starts  to  re-freeze.  There  is  a gap  between 
the  melt  and  freeze  fronts,  signifying  the  storage  of  the  latent 

the  cycle  when  the  heat  is  withdrawn  from  the  medium.  This 
retracted  front  eventually  merges  with  the  advancing  re-freeze 
front,  resulting  in  the  disappearance  of  the  latent  heat  stored. 

heat  at  the  end  of  the  cycle.  However,  for  the  medium  imposed 


problems  in  cylindrical  i 


Kolodner  [16]. 


4.  It  is  noted  that  a constant  element  has  been  used  in  the 
development  of  the  SS1I  and  BEM  in  this  work.  It  is  desirable  to 


SAMPLE  SSM  FORTRAN  PROGRAM  FOR  CYCLIC 
TEMPERATURE  BOUNDARY  CONDITION  WITH  SUBCOOLING  OR  SUPERHEATING 


VRITE(16,*) ’ INPUT  [ 
VRITE(16,122)TM 
VRITE(16,133)CK 
VRITE(16,144)RH0 
VRITE( 16, 155)CP 
VRITE(16,166)CL 
VRITE( 16 , 177)DELTM 
WRITE(16,188)TRTH 
VRITE(16,  “ 


VR1TE(16," 

WRITE(16,. 

WRITE(17,< 


INTERNAL  COMPUTATION  S' 


RX(l)=O.ODO 
TIN( I)=O.ODO 
CONTINUE 


E INTERFACE  POSITION  INTERVAL  Rl(NF) 


VRITE(.,.)TFF,’ 

RA=R1(NF-I)/2.0I 


R1(NF)=ZBRENT(TFUNC1  ,RA,R8,T0L) 
WRITE(.,.)  NF.Rl(NF) 
VRITE(16,1000)TFF,R1(NP) 

—UPDATE  INTERFACE  POSITION  VELOCITY 


RV1(NF)=(R1(NP)-R1(NF-1))/DELTM 


— COMPUTATION  OF  TEMPERATURES  AT  INTERIOR  PI 


RX(1)=RIX 

TIN(1)=FX(TIHE) 

VRITE(17, 199)TIME.RIX,TIN 

RIX=R1(NF)/10. ODO.I 

RX(U1)=RIX 

CALL  INTERNAL(TINTL) 

CONTINUE 


RIX=R1(NF)+R1(NF)*12. ODO.I 
E IF((NF  -GT.  3)  .AND.  (NF  .1 
RIX=R1 (NF)+R1 (NF)*5. ODO.I 
IE  IF((NF  .GT.  6)  .AND.  (NF  .1 
RIX=Rl(NF)-fRl(NF).3. ODO.I 
IE  IF( (NF  -GT.  12)  .AND.  (NF  . 
RIX=R1(NF)+RI(NF).2. ODO.I 
IE  IF((NF  .GT.  20)  .AND.  (NF  . 
RIX=R1(NF)+R1(NF).I 
IE  IF((NF  .GT.  20)  .AND.  (NF 
R1X.R1 (NF)tRl (NF)/2. ODO.I 

RIX=Rl(NF)+Rl(NF)/3. ODO.I 


L INTERNAL(TINTL) 


-COMPUTATION  OF  STORED 


-COMPUTE  TVO  INTERFACE  POSITIONS  R1.R2 
-INITIAL  GUESS  OF  INTERFACE  POSITION  R2(NF) 


SUCCESSIVE  ITERATION  FOR 


VRITE(.,.)NF,’  ENTER  ASSUMED  INITIAL  R2=?' 
READ(»,.)R2(NF) 

R20LD=R2(NF) 

R2(NF)=R2(NF-1 ) 


Rl(NF)=ZBRENT(TFUNCI,RA,RB,TOL) 

R10LD=R1(NF) 

-SOLVE  FOR  R2(NF) 


RB=5.0D0.R2(NF-1) 

R2(NF)=ZBRENT(TFUNC2,RA,RB,T0L) 

R1(NF)=ZBRENT(TFUNC1  ,RA  ,RB,TOL) 

-CHECK  ACCURACY  OF  COMPUTED  R1,R2 

FLAG1=DABS((R1(NF)-R10LD)/RI(NF)) 

FLAG2=DABS((R2(NF)-R20LD)/R2(NF)) 

IF( (FLAG  1 .LT.  EPS)  .AND.  (FLAG2  .LT.  EPS))  THEN 


VRITE(.,.jNF,Rl(NP),R2(NF) 
VRITE(  16 , 1200)TFF,R1(NP)  ,R2(NF) 


--UPDATE  INTERFACE  POSITION  VELOCITY 

RV1(NF)=(R1(NF)-R1(NP-1))/DELTM 
IF(NF  . EQ . MT+1)  THEN 
RV2(NF)=R2(NF)/DELTM 

^RV2(NF)=(R2(NF)-R2(NF-1))/DELTM 
—COMPUTATION  OF  TEMPERATURES  A 
IF(MOD(NF,NPS)  .NE.  0)  GO  TO 


INTERIOR  POINTS 


RIX=R2(NF)/10.0D0*I 
CALL  INTERNAL(TINTL) 

CONTINUE 
VRITE(1T,*)’  ’ 

DO  310  1=1,10 

RIX=(R1  (NF)  -R2(NF)  )/10. 01 
RX(I+11)=RIX 
CALL  INTERNAL(TINTL) 
TIN(I*11)=TINTL 
CONTINUE 

RIX=R1(NF)*R1(NF)/2.0D0.1 
CALL  INTERNAL(TINTL) 

TIN( 1+21 )=TINTL 
CONTINUE 


VR1TE( • , ' ITERATION  STEP=’,NSS 

WRITE( • ,•) 'R1 ' ,R1(NF) .RIOLD 
VRITE(. , . ) ’ R2 ' ,R2(NP) ,R20LD 

R10LD=R1(NF) 

R20LD=R2(NF) 


IF(MPP  .EQ.  0)  GO  TO  100 
-COMPUTE  INTERIOR  TEMPERTURE  OF  MERGED  PHASE 
IF(MOD(NF,NPS)  -NE.  0)  GO  TO  100 

VRITE(17*»)* MERGE ’ 

TIME=TM1<NF*DELTM 

TIN(1)=FX(TIME) 

VRITE(17,199)TIME,RIX,TIN(1) 

RIX=R2(MPNF)/3.0D0.I 
CAUL  ONEPHASE(TINTL.MPNF) 

CONTINUE 


CALL  ONEIIEATS 
100  CONTINUE 


FORMAT( // ’ MELTI NG  TEMPERATURE’ ,1 
FORMAT( ’THERMAL  CONDUCTIVITY^ ,1 
FORMAT( ’DENSITY=’ ,F7.2) 

FORMATf ’SPECIFIC  HEAT  CAPACITY= 
FORMAT( ’ LATENT  HEAT= ’ , F10 . 2 ) 
FORMAT(/’TIME  STEP  SIZE=’.F5.3> 
FORMAT( ’TRANSI 
F0RMAT(2X,F:" 

0 F0RMAT(3X,F: 

0 F0RMAT(3X,P: 


3//) 


SUBPROGRAM  FUNCTION 


FUNCTION  TFUNCl(X) 

IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 
PARAMETER(  MAXF=2500  ,PI=3. 11 1592653589793D0  ) 


IF(NF  .EQ.  1)  THEN 

RV1(NF)=R1(NF)/(TFF2-TFFJ) 

ELSE 

RV1(NF)=(R1(NF)-R1(NF-1))/(TFF2-TFP1) 


IF(THI  .Eq.  O.ODO)  THEN 
TF2=TMI 

CALL  RDMBERG(FVNC1,TF1,TP2,QINT) 


CALL  R0MRERG(FUNC1 ,TF1 ,TF2,A) 


TF1=TMI+(1NF-1)*DELTM 
AA=DSQRT(TFF2-TF1 ) 
BB=DSQRT(TFF2-TF2) 

CALL  R0MBERG(FUNC2A,AA,I 
CALL  ROMBERG(FUNC2B,AA,l 
B1=RV1(INF)*B1 
B2=RV1(INF).B2 


IF(NF  .LE.  HT)  THEN 
SUMC=O.ODO 

RV2(NF)=R2(NF)/(TFF2-TFF1 ) 


RV2(NF)=(R2(NF)-R2(NF-1))/(TFF2-TFF1) 


S.qiNT) 


RV1  (NF)=(R1(NF)-R1(NF-1 ) )/ (TPF2-TPF1 ) 


TfurC2=(R2(NF)*(SUMA+QINT)-CL/CP*SUMB+CL/CP*SUIIC) 
S /2.0!I0/DSQRT(PI-ALPBA)-TR 


PARAMETER(PI=3. 141592653589793D0) 
FX=20 . 090.20. ODO.DCOS(PI/12-O.X.PI ) 


X).DEXP(-RICNF)..: 

>QRT(TFF-X)»«3 


(NF)..2/4 .ODO/ALPHA/(TFF-X) ) 


X /T/NF.INF.TMI.TMF.DELTSt 

N /C/ALPHA 

N /R1/R1(MAXF),RV1(MAXF) 


TP1=T*H(  INF-1  ).DELTH 
IF(INF  -EQ.  1)  THEN 
A=RV1(INF) 

B=(R1(!nF-  1 ).TF2-R1  ( INFJ.TFl  )/(TF2-TPl ) 
RX=A.(TPF-X..2)*H 
IFFUNraA'0°0D0>) 

ELFUNC2A=-2.0D0*DEXP(-(Rl(NF)tRX)«*2/4.0D0/ALPHA/X/X) 


CQHMQN^/T/NF^NF^TNI  ,THF,DELTU 

COMMON  /C/ALPIIA 

COMMON  /R1/R1(MAXF),RV1(MAXF) 


TF2=TM1+INF»DELTM 

TFl=TMI+(  INF-1J-DELTM 
IF(1NF  .EQ.  1)  THEN 
A=RV1(INF) 

B=(R1(  INF-1  )*TF2-R1(INF)*TP1)/(TF2-TF1) 
RX=A.(TFF-X*.2)+B 

IF(X  .EQ.  O.ODO)  THEN 
FUNC2B=-2 . ODO 

ELFUNC2B=-2.0D0.DEXP(-(R1(NF)-RX)..2/4.0D0/AEPHA/X/X) 

RETURN 


IN  FUNC3A(X) 

:T  DOUBLE  PRECISION 


B=-R2( INFJ.TFl /(TF2-TF1 ) 

B=(R2(InF-1).TF2-R2(INF).TF1)/(TF2-TI 


P(-(Rl(NF)+RX).-2/4.0D0/ALPHA/X/X) 


FUNCTION  FUNC3B(X) 

IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 
PARAMETER(MAXF=2500) 

COMMON  /T/NF,INF.TMI.TMF,DELTM 

COMMON  /C/ ALPHA 

COMMON  /R1/R1(MAXF),RV1(MAXF) 


COMMON  /R2/R2(MAXF),RV2(MAXF),MF,I 


1F(INF  .Eq.  MT*1)  THEN 
A=RV2(INF) 

B=-R2(  INF).TF1/(TF2-TF1 ) 

IF(X  .Eq.  O.ODO)  THEN 
IF(R1(NF)  .Eq.  R2(NF))  THEN 
FUNC3B=-2.0D0 

FUNC3B=0.0D0 
END  IF 

BLF«NC3B=-2.0DO.DEXP(-(R1(NF)-RX)«.2/4.0DO/ALPHA/X/X) 


T DOUBLE^PRECISIOI 


IN  /T/NF,1NF,TMI,TMF,DELTM 
IN  /C/ALPHA 

IN  /R2/R2(MAXF),RV2(MAXF) ,M! 


IF((X-TFF)qjEq.  O.OI 


FUNCTION  FUNCSA(X) 

IMPLICIT  DOUBLE  PRECISION(A-H,0- 
PARAMETER(MAXF=2500) 

COMMON  /T/NF, 1NF.TMI .TMF.DELTM 
COMMON  /C/ALPHA 

COMMON  /R2/R2(MAXF),RV2(MAXF),MI 


A=RV2(INF) 

B=-R2(INF)*TF1/(TF2-TF1) 

A=RV2(INF) 

B=(R2( INF- 1 ).TF2-R2( INF).TF1 )/(TF2-TFl ) 
RX=A*(TPF-X..2)tB 


'( - (R2(NF)+RX)*»2/ 4 .ODO/ ALPHA /X/X) 


^precision 


IF( INF  .EQ.  MT+1)  THEN 
B=-R2(1NF).TF1/(TF2-TF1) 


B=(R2(  INF-1  ).TF2-R2(  INF).TF1  )/(TP2-TF: 
RX=A»(TFF-X**2)+B 


IF(X  .EQ.  O.ODO)  THEN 
FUNC5B=-2.0D0 


(-(R2(NF)-RX)« 


2/4 . ODO/ ALPHA /X/X) 


FUNCTION  FUNC6A(X) 

IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 
PARAMETER( MAXF=2500 ) 

COMMON  /T/NF,1NF,TMI .TMF.DELTM 

COMMON  /C/ALPIIA 

COMMON  /R1/RI(MAXF),RV1(MAXF) 

COMMON  /R2/R2(MAXF) ,RV2(MAXF) ,MF,MT 


IF(INF  .EQ.  1)  THEN 
A=RVJ(INF) 

BrO.ODO 

A=RV1(INF) 

B=(R1(INF-1).TF2-R1(1NF).TF1)/(TF2-TF1) 


ELFUNC6A=-2.0DO.DEXP(-(R2(NF)*RX)..2/4.0D0/ALPHA/X/X) 


FUNCTION  FUNC6B(X) 

IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 
PARAMETER(MAXF=2500) 

COMMON  /C/ ALPHA 

COMMON  /R1/R1(MAXF),RV1(*AXP) 

COMMON  /R2/R2(MAXF),RV2(MAXF),MF,MT 


TF1=TMI.(INF-1).DELTM 

A=RVI(INF) 

A=RV1(INF) 

B=(R1 ( INF-1 J.TF2-R1 ( INFJ.TFl )/(TF2-TFl ) 
RX=A.(TFF-X..2)tB 

IF(R1(NF)  ?Eq!  R2(NF))  THEN 
FUNC6B=-2.0D0 

FUNC6B=0.0D0 


FUNC6B=-2.0D0*DEXP(-(R2(NF)-RX)*«2/4.0D0/ALPHA/X/X) 

RETURN 


SUBPROCRAM 

USING  MODIFIED  BISECTION  METHOD  (BRENT'S  METHOD) 
TO  FIND  THE  ROOT  OF  FUNCTION  TFUNC 


FUNCTION  ZBRENT(TFUNC,Xl,X2,TOL) 


gs 

If— 

"Ss&'m'rm  -*nd-  d 

"ESSE* 


SUBPROGRAM 
ROMBERG  INTEGRATION 

INPUT— A, B, EPS  (INTERVAL  AND  CRITERION) 


T(l,l)=(B-A).(FUNC(A)*FUNC(B))/2,ODO 

T(I,2)sT(I,l)/2.ODO+(B-A)«FUNC((A+B)/2.OD0)/2.Ol 

T(2,1)=(4.0D0.T(1,2)-T(l,l))/3.0D0 


—SUCCESSIVE  APPLICATION  OP  TRAPEZOIDAL  Rl 

DELX=(B-A)/2,0D0..(J-1) 

X=A-DELX 

N=2..(J-2) 


SUM=SUM+FUNC(X) 

I CONTINUE 

T(l» J)=T(I , J-1)/2.0IX 
— EXTRAPOLATION 


T(L,K)=(4.0D0..(L-1).T(L-1, 
* (I.ODO..(L-1)-1.0DO; 

I CONTINUE 

— CBECK  ACCURACY  CRITERION 
IF(T(J,1)  .EQ.  O.ODO)  THEN 


FUNCTION  FINL2A(X) 

IMPLICIT  DOUBLE  PRECISION(A-n.O-Z) 
PARAMETER(MAXF=2500) 

COMMON  /T/NF,INF,TMI,HIF,DELTM 
COMMON  /C/ALPUA 
COMMON  /R1/R1(MAXF).RV1(MAXF) 
COMMON  /RINT/RIX 


TF2=TMI+I  NF*DELTJJELT1( 

IF(INF  .EQ.  1)  THEN 
ArRVl(INF) 

ELSE 

A=RV1(INP) 

B=(R1 ( INF-I J.TF2-R1 ( INF)«TF1 )/ (TF2-TF1 ) 
END  IF  x 2)  B 


IF(X  .Eq.  O.ODO)  T 


2 . ODO-DEXP( - ( R1 X.RX ) ..2/4 . ODO/ALPHA/X/X) 


IF(X  .EQ.  O.ODO)  THEN 

IF(RIX  .EQ.  Rl(NF))  THEN 

ELSE 


( - (RIX-RXJ.-2/4 .ODO/ALPHA/X/X) 


RX=A.(TFF-X.*2)*B 
IF(X  .EQ.  O.ODO)  THEN 
IF(RNL3BEQ2  OHO1"1’’ 
FINL3B=O.ODO 


(-(RIX-RX)**2/4.0D0/ALPHA/X/X) 


CALL  R0MHERG(FINL2B,A. 
B1=RV1(INF).B1 
B2=RV1(INF).B2 
SUHB1=SU1'"'  ” 


3RG(FINL3B,AA,BB,C2) 
C1=RV2(INF).C1 
C2=RV2(1NF).C2 
SUMC1=SUMC1*C1 


10, CP, CL, TM 

IN  /Rl/Rl ( MAXF ) , RV1 ( MAXF ) 

IN  /R2/R2(HAXF) ,RV2(RAXF) ,HI 

COMMON  /HEAT/RX(51 ) ,TIN(51) 
--COMPUTE  THE  STORED  LATENT  HEAT 
SLH=RH0.(R1(NF)-R2(NF) ).CL 


• ^ODO*TIM21  j+4  0DO"TIN(22)*2.0DO"TIN(23)*4.0DO-TIN(2'l)* 
S hJImoSiN  SU2.0DO.TIN(27)*4.0DO.TIN(28  ♦ 

S 2.ODO.TIN(29)*4.0DO.TIN(30)*2.ODO.TIN(31)»4.ODO.TIN(32)+ 

S 2*0DO*TIN(33)+4.Or  * " 

S TIN(41)j 


SUBPROGRAM 
MERGING  START 
STORED  BEAT  COMPUTATION 

SUBROUTINE  ONEHEATS 
IMPLICIT  DOUBLE  PRECISION(A- 
PARAMETER(MAXF=I 
COMMON  /T/NF,INP,l»i.i»r, UN- 
COMMON /Cl/CK,RBO,CP ,CL,TM 
COMMON  /R1/RI(MAXF),RV1(MAXF) 


1=3. 14I592653589793D0) 


COMMON  /R2/R2(MAXF),RV2(MAXF),MF,MT 
COMMON  /HEAT/RX(51),TIN(51) 

COMPUTE  TUB  STORED  SENSIBLE  BEAT  IN  THE  MATERIAL 
(SIMPSON’S  1/3  RULE) 


TOTAL=SSH 

VRITE(  17 ,200)TOTAL , SSH 
F0RMAT(/1X, ’STORED  HEAT=\E14.7,2X,2X 
S, ’SENSIBLE:’ ,E14.7) 


APPENDIX  B 

SAMPLE  BEM  FORTRAN  PROGRAM  FOR  CONSTANT 
OR  MONOTONIC  TEMPERATURE  OR  FLUX  CONDITIONS  OF  TIME 


VRITE(16,.) 'THERMAL  CONDUCTIVITY  OF  PHASE  1 =’,CK1 
VRITE( 16,.) 'THERMAL  CONDUCTIVITY  OF  PHASE  2 =\CK2 
VRITE(16,.) ’DENSITY.’, RHO 

VRITEf  16,.) ’SPECIFIC  HEAT  CAPACITY  OF  PHASE  1 =’,CP1 
VRITE(  16, . ) ’SPECIFIC  HEAT  CAPACITY  OF  PHASE  2 = ’,CP2 
VRITE( 16,.) 'LATENT  HEAT.’, CL 

VRITE(16,.)’ 

VRITE(16,.)'  TIME  INTERFACE’ 

ALPHAI.CK1/RHO/CP1 


ALPHA2.CK2/RHO/CP2 
MF=INT((TMF-TMI)/DELTM) 
CALL  INPUT(MF ,TMI .DELTM) 


R(NF)=ZBRENT(TFUNC,Rl,R2,TOL) 
CALL  OUTPUT 


100  CONTINUE 


COMMON  /R/R(MAXF) ,RV(MAXF) 
COMMON  /S/H(2,2),G(2,2),DBC(2,MI 
COMMON  /B/K0DE(2,MAXF),BC(2,MAXi 
EXTERNAL  FTMP 1 , FTMP2 , FtJMP 
R(NF)=X 


IF(NF  .Eq.  1)  THEN 

RV(NF)=R(NF)/(TFF2-TFF1 ) 

RV(NF)=(R(NF)-R(NF-1))/(TFF2-TPF1) 

END  IF 

COMPUTE  THE  PLUX  OF  INTERFACE  OF  PHASE  1 

S0L1=(0(2.2)-DBC(1,NF)-C(1,2).DBC(2,NF))/ 
S <G(2,2).G(I,1)-0(1.2).0(2,1)) 

SOL2=(G(2,l).DBC(l,NF)-G(l,l)-DBC(2,NF))/ 
$ (G(2,l)-G(l,2)-G(l,l)-G(2,2)) 

DBC(l,NF)=SOLl 
DBC(2,NF)=SOL2 

--COMPUTE  THE  FLUX  OF  INTERFACE  OF  PHASE  2 
SUMA=O.ODO 

CALL  I NTE22(  ALPHA2 , HSUM , GSUM , I NF) 

IF(INF  -EQ.  NF)^THEN 

SUMA=SUMA+(HSUM-RV(INF)*GSUM/ALPHA2)»TM 


I=DSQRT(ALPBA2/Pi; 


I.0D0/CK2 


TFUNC=CK2*QM(NF)-CK1*DBC(2,NF)-R1 


IMPLICIT 


« FTMP2(X,DUM1,IDUH) 

T DOUBLE  PRECIS  I ON(A-B.O-Z) 

ER(MAXF=500 , PI  =3 . 14 1592653589793D0 ) 

COMMON  /T/NF ,TMI  ,TMF , DELTM 
COMMON  /C/ALPHA  1.ALPHA2 
COMMON  /R/R(MAXF),RV(MAXF) 

TFF=TMI*NF.DELTM 

(TMi-X)/DSQRT(4.0b0.AI 

•ERF(ARG2))/(DSQRT(TFF- 


N /T/NF,TMI,TMF,DELTM 
OMMON  /C/ALPBA1 ,ALPBA2 
OMMON  /R/R(MAXF),RV(MAXF) 

RGl=R(NF)..2/4-0D0/ALPBA2/(TFF-X) 

RG2=R(NF)»DSQRT(TMI-X)/DSQRT(4.0D0*ALPBA2*(TFF-TM. 


X)) 

|MP=BCG(X).DI 


K)+ERF(ARG2))/DSQRT(TFF-X) 


EQUATION 


„u.l(A-H,0-Z) 

00, Pl=3 . 141592653589793D0) 
;SI0N'VH(2,NAXF),V0(2,MA*F),VVB(2),VVG(2) 

IN  /T/NF ,TMI  ,TMF,  DELHI 
IN  /C/ALPBA1 , ALPBA2 
IN  /B/K0DE(2,MAXF) ,BC(2,MAXF) 


—INITIALIZE  VVB.VV 
VVR(I)=0.0D0 
— TIME  ITERATION  FI 


,L  INTEGRATION  SUBROUTINES 

((I  .EQ.  1)  .AND.  (J  .EQ.  1))  TREN 

CALL  INTE1I(ALPBA1,B(I,J),G(I,J),INF) 

IF  (NF  .EQ.  INF)  TflEN 
B(I,J)=B(I,J)-0.5D0 

5E^' IF* ( ( I -EQ.  1)  .AND.  (J  .EQ.  2))  TBEN 

CALL  INTEI2(ALPBA1,B(I,J),G(I,J),INF) 

SE  IF  ((I  .EQ.  2)  .AND.  (J  -EQ.  1))  TBEN 

CALL  INTE21 ( ALPHA1 , B( I , J) ,G( 1 , J ) , INF) 


. CONTINUE 
—COMPUTE  B,( 


H(I,1)=B(I,1) 

B ( 1 ,2)=RV(INF)/ALPBA1»G(1 ,2)-H(I,2) 

G(I,1)=G(I,1) 

600  G(I,2)=-G(I,2) 


IF(INF  .EQ.  NP) 


I CONTINUE 

—COMPUTE  THE  EAST  TIME  STEP’S 

10  VVH(I)=VVH(I)-VVG(I) 

— REORDER  THE  SYSTEM  OF  EQNS  T( 

IF(KODE(J,NF))1100, 1100,1200 

CH=G(1 ,J)~  ' 

G(I,J)=-H(I,J) 

•0  H ( I ,J)=-CH 
10  CONTINUE 

—COMPUTE  RHS  KNOVN  VECTOR  DBC 

Sar 


COMMON  /R/R(MAXF) ,RV(MAXF) 


I ) -DSQRT(TFF-TF2) ) 


issrsasB.,., 

PARAMETER) MAXF=500 ,P1=3. 141 5926535881 


'lsTMl+(INP-l)«DELTH 


2-R(INF)*TFl  )/(TF2-TF: 


G=AL/2.0DO/A*(1.0D0-ERF(Cl)+2.0D0*H) 

A2=(2.0DO»A»TFF+B-A*TF2)/DSQRT(4.0DO*AL« 
AU(2.0D0.A.TFF*B-A.TF1)/DSQRT(4.0D0.AL< 
CALL  RDHBERC(PH12,A1,A2,H,AL,INF) 
C2=(A.TF2*B)/DSQRT(4.0D0.AL.(TFF-TF2)) 
C1=(A.TF1+B)/DSQRT(4-0D0.AL.(TFF-TF1)) 
CALL  R0MBERG(FG12,CI,C2,0,AL,INF) 


FUNCTION  FH12(X,AL,INF) 

IMPLICIT  BOUBLE  PRECISION(A-H,0-Z) 
PARAMETER(MAXF=500 , PI =3 - 14 1592653589793D0 ) 
COMMON  /T/NF,THI,TMF,DELTM 
COMMON  /R/R(MAXF),RV(MAXF) 


A=RV(INF) 

B=-R(INF).TF1/(TF2-TF1) 

A=RV(INF) 


:(R(INF-1 ).TF2-R( INFJ.TFl )/(TF2-TFl ) 

I .Eq.  O.ODO)  THEN 
I12=-DEXP( -X..2)/DSQRT(PI ) 

I12=-DEXP(-X..2*A.(A.TFP*B)/AL)/BSQRT(PI) 


FG12=AL*DEXP(-X**2)/A/DSQRT(PI ) 


1F(NF  .Eq.  INF)  THEN 

AUR(NF)/BSQRT(4.0D0.AL.(TFF-TF1)) 

G=R(NF)/2.0B0/DSqRT(PI).(DEXP(-A1..2)/Al-DSQRT(PI). 


giggSgi 

■ Ssirff 

^DR=(ERF(A2)-ERF(Ali 


0D0.AL.(TFF-TF2)) 


[=3.141 5926535897 93D0) 


IF(INF  .EQ.  1)  THEN 
A=RV(INF) 

B=-R(  INP)»TP1  /(TF2-TP1 ) 
A=RV(INP) 


1F(A  -El).  O.ODO)  THEN 

G=DSQRT(AL/PI ).DSQRT(TFF-TF1) 

A2=(R(NF)-R( INF) )/DSQRT(4 . OBO.AL.(TFP-TF2) ) 
Al=(R(NF)-R(INF))/BSqRT(4.0D0.AU(TFF-TFl)) 

CALL  ROMBERG (FG22,A1 ,A2>G,AL,INF) 
G=(R(NP)-R(INF))/2.0D0/DSQRT(PI) 

S .(DEXP(-A1..2)/A1-DEXP(-A2..2)/A2-2.0DO.G) 

CALL  R0MBERG(FH22,A1 ,A2,H,AL,INF) 

IF(NF  .EQ.  INF)  THEN 

B1=(R(NF)-2»A»TFF+A»TF1-B)/DSQRT(4.0D0*AL»(TFF-TF1)) 

C1=(R(NF)-A.TF1-B)/BSQRT(4.0D0.AL.(TFF-TFI)) 

G=AL/2 .ODO/A. ( ERF(C1 )*2 . ODO.H) 

B2=(R(NF)-2.A.TFF*A.TF2-B)/DSQRT(4.0D0.AL.(TFF-TF2)) 
B1=(R(NF)-2.A.TFF*A.TF1-B)/DSQRT(4.0D0.AL.(TFF-TF1» 
CALL  ROMDERG ( FH22 , HI , B2 , II , AL , INF) 
C2=(R(NF)-A.TF2-B)/DSQRT(4.0DO.AL.(TFF-TF2)) 
C1=(R(NF)-A«TF1-B)/DSQRT(4.0DO.AL.(TFF-TF1)) 

CALL  R0MBERG(FG22,CI,C2,G,AL,INF) 


FUNCTION  FH22(X,AL,INF) 

IMPLICIT  DOUBLE  FRECISION(A-H.O-Z) 
PARAMETER MAXF=500 , P I =3 . 14 1592653589793D0 ) 
COMMON  /R/R(MAXF),RV(MAXF) 


IF(INF  .EQ.  1)  THEN 
B=-R( INF).TF1/(TF2-TF 


TF2-R( INK).TF1 )/(TF2-TI 

K))  THEN 
;—2)/DSqRT(PI) 
:..2-A.(R(NF)-A. 


P-B)/AL)/DSQRT(PI) 


FUNCTION  FG22(X,AL,INF) 

IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 
PARAMETER(MAXF=500,PI=3. H1592653589793D0) 

COMMON  /R/R(MAXF),RV(MAXF) 


IF(INF  -El).  1)  THEN 
A=RV(INF) 

B=-R(INF).TF1/(TF2-TI 


F2-R(INF).TF1)/(TF2-TI 


■L*DEXP(-X»»2)/A/DSQRT(PI ) 


SUBPROGRAM 

ROMBERG  INTEGRATION  PROGRAM 

INPUT— A, B, EPS  (INTERVAL  AND  CRITERION) 


DELX=(B-A)/2.0D0..(J-1) 

N=2..(J-2) 

o JSS.« 

"SuLMiSir”"”  ” 

IFjDJK((I(J,l)-T(J-l,l))/T(J,l»  .GE.  EPS)  THI 


FUNCTION  ZBIiENT(TFUNC,Xl,X2,T0L) 


IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 
EXTERNAL  TFUNC 

PARAMETER  (ITMAX=100,EPS=3.D-16) 


O.DO)  PAUSE  'Root  oust 


IF(OABS(FC).LT.DABS(FB»  Tl 


■AND.  DABS(FA).GT.DABS(FB))  THEN 


R)-(B-A).(R-l.DO)) 


IF(P.GT.O.DO)  Q=-Q 
P=DABS(P) 

IF(2.D0.P  -LT.  DMIN1 (3 . 


•Q-DABS(TOLI.q),DABS(E.q>)) 


B=B*DSIGN(T0L1  ,XM) 

FBrTFUNC(B) 

CONTINUE 


FUNCTION  ERF(X) 

IMPLICIT  DOUBLE  PRECISION(A-H, 
IF(X.LT.O.ODO)TBEN 

ERF=-GAMMP(-5D0,X..2) 


0. DO) (PAUSE  ’NEGATIVE  A 


FUNCTION  GAMMQ(A,X) 

IMPLICIT  DOUBLE  PRECISION(A-B.O-Z) 
IF((X.LT.O.ODO).OR.(A.LE.O.ODO))PAUSE  ’NEGATIVE  A 
IF(X.LT.A+l.ODO)THEN 
CALL  GSER( GAMSER , A , X , GLN) 

GAMMQ=1 . DO-GAMSER 

CALL  GCF(GAMMCP,A,X,GLN) 


SUBROUTINE  GCF(GAMMCF,A,X,GLN) 
IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 
PARAMETER  (ITMAX=l00,EPS=3.D-9) 
GLN=GAMMLN(A) 


IF(DABS((G-GOLD)/G).LT.EPS)GO  Tl 

GOLDsG 

CONTINUE 


iLSVRITE(16,350)l,BC(I,NF),DBC(I,NF),«M(NF) 


F0RM»T(«,’1 


300  FORMAT!  M.2> 
350  FORMAT! 14, 2> 
400  FORMATflX.F! 


/•ATSiS.'teSfSSfcE  '■■ 
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